## Demonstration of using R for analyzing data as in the paper:
# You need: 
# R v. 2.13.0 or above (tested only under Windows)
# library nlme (see http://cran.r-project.org)
# Your data set 

# Read your file and exploration of raw data:
dat <- read.csv ("Plividus.csv", sep=";", header=T)
str (dat)
summary (dat)
plot (log10 (dat$gdw) ~ log10 (dat$diameter))  #where gdw is the gonad dry weight
hist (dat$dgw)
qqnorm (dat$dgw)

# linear regression 1: 
model1<- lm (log10 (gdw) ~ log10 (diameter), data=dat)
# check residuals:
par (mfrow =c(2,2))
plot (model1)

# linear regression 2:
model2<- lm (log10 (gdw) ~ log10 (diameter-27.9), data=dat)
# check residuals:
par (mfrow=c(2,2))
plot (model2)

# We compare both models by AIC:
AIC (model1, model2)

# We define the SGI and the LGI:
dat$sgi <- residuals (model1) 
dat$lgi <- dat$gww /dat$iww *100  # where gww = gonad wet weight and iww: Individual wet weight

# We analyze the effect of the diameter on the SGI and LGI using linear mixed models:
library (nlme) 

mixmod1<-lme (sgi ~ diameter, random=~1|month, corr = corAR1(form=~1), weights = varIdent (form = ~1|month), data=dat)
summary (mixmod1)

#Check residuals: 
a<- resid (mixmod1, type="normalized")
b<- fitted (mixmod1)
plot (b, a, xlab="Fitted values", ylab="Residuals")
plot (dat$diameter, a, xlab="Diameter", ylab="Residuals")
plot (dat$month, a, xlab="Month", ylab="Residuals")
qqnorm (a);qqline (a)
hist (a)
plot (ACF (mixmod1, resType = "normalized"), alpha=0.05)
plot (augPred (mixmod1, level = 0:1, primary=~diameter), layout = c(5,5), cex=0.5, main="", ylab="SGI", xlab="Diameter")     

mixmod2<- lme (lgi ~ diameter, random=~1|month, corr = corAR1(form=~1), weights = varIdent (form=~1|month), data=dat)
summary (mixmod2)

#Check residuals: 
a<- resid (mixmod2, type="normalized")
b<- fitted (mixmod2)
plot (b, a, xlab="Fitted values", ylab="Residuals")
plot (dat$diameter, a, xlab="Diameter", ylab="Residuals")
plot (dat$month, a, xlab="Month", ylab="Residuals")
qqnorm (a); qqline (a)
hist (a)
plot (ACF (mixmod2, resType = "normalized"), alpha=0.05)
plot (augPred (mixmod2, level = 0:1, primary=~diameter), layout = c(5,5), cex=0.5, main="", ylab="LGI",  xlab="Diameter")